namespace Eigen {

namespace internal {

    template <typename FunctorType, typename Scalar>
    DenseIndex fdjac1(const FunctorType& Functor,
                      Matrix<Scalar, Dynamic, 1>& x,
                      Matrix<Scalar, Dynamic, 1>& fvec,
                      Matrix<Scalar, Dynamic, Dynamic>& fjac,
                      DenseIndex ml,
                      DenseIndex mu,
                      Scalar epsfcn)
    {
        using std::abs;
        using std::sqrt;

        typedef DenseIndex Index;

        /* Local variables */
        Scalar h;
        Index j, k;
        Scalar eps, temp;
        Index msum;
        int iflag;
        Index start, length;

        /* Function Body */
        const Scalar epsmch = NumTraits<Scalar>::epsilon();
        const Index n = x.size();
        eigen_assert(fvec.size() == n);
        Matrix<Scalar, Dynamic, 1> wa1(n);
        Matrix<Scalar, Dynamic, 1> wa2(n);

        eps = sqrt((std::max)(epsfcn, epsmch));
        msum = ml + mu + 1;
        if (msum >= n)
        {
            /* computation of dense approximate jacobian. */
            for (j = 0; j < n; ++j)
            {
                temp = x[j];
                h = eps * abs(temp);
                if (h == 0.)
                    h = eps;
                x[j] = temp + h;
                iflag = Functor(x, wa1);
                if (iflag < 0)
                    return iflag;
                x[j] = temp;
                fjac.col(j) = (wa1 - fvec) / h;
            }
        }
        else
        {
            /* computation of banded approximate jacobian. */
            for (k = 0; k < msum; ++k)
            {
                for (j = k; (msum < 0) ? (j > n) : (j < n); j += msum)
                {
                    wa2[j] = x[j];
                    h = eps * abs(wa2[j]);
                    if (h == 0.)
                        h = eps;
                    x[j] = wa2[j] + h;
                }
                iflag = Functor(x, wa1);
                if (iflag < 0)
                    return iflag;
                for (j = k; (msum < 0) ? (j > n) : (j < n); j += msum)
                {
                    x[j] = wa2[j];
                    h = eps * abs(wa2[j]);
                    if (h == 0.)
                        h = eps;
                    fjac.col(j).setZero();
                    start = std::max<Index>(0, j - mu);
                    length = (std::min)(n - 1, j + ml) - start + 1;
                    fjac.col(j).segment(start, length) = (wa1.segment(start, length) - fvec.segment(start, length)) / h;
                }
            }
        }
        return 0;
    }

}  // end namespace internal

}  // end namespace Eigen
